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ABSTRACT 



Context. Simple models fail to describe the observed spectra of X-ray dim isolated neutron stars (XDINSs). Interpretation of these 
spectra requires detailed studies of radiative properties in the outermost layers of neutron stars with strong magnetic fields. Previous 
studies have shown that the strongly magnetized plasma in the outer envelopes of a neutron star may exhibit a phase transition to 
a condensed form. In this case thermal radiation can emerge directly from the metallic surface without going through a gaseous 
atmosphere or, alternatively, it may pass through a "thin" atmosphere above the surface. The multitude of theoretical possibilities 
complicates modeling the spectra and makes it desirable to have analytic formulae for constructing samples of models without going 
(— I \ through computationally expensive, detailed calculations. 

Aims. The goal of this work is to develop a simple analytic description of the emission properties (spectrum and polarization) of the 
condensed, strongly magnetized surface of neutron stars. 

Methods. We have improved the method of van Adelsberg et al. (2005) for calculating the spectral properties of condensed mag- 
netized surfaces. Using the improved method, we calculate the reflectivity of an iron surface at magnetic field strengths B ~ 10 12 G 
- 10 14 G, with various inclinations of the magnetic field lines and radiation beam with respect to the surface and each other. We 
construct analytic expressions for the emissivity of this surface as functions of the photon energy, magnetic field strength, and the 
three angles that determine the geometry of the local problem. Using these expressions, we calculate X-ray spectra for neutron stars 
with condensed iron surfaces covered by thin partially ionized hydrogen atmospheres. 

Results. We develop simple analytic descriptions of the intensity and polarization of radiation emitted or reflected by condensed iron 
surfaces of neutron stars with strong magnetic fields typical for isolated neutron stars. This description provides boundary conditions 
at the bottom of a thin atmosphere, which are more accurate than previously used approximations. The spectra calculated with this 
improvement show absorption features different from those in simplified models. 

Conclusions. The approach developed in this paper yields results that can facilitate modeling and interpretation of the X-ray spectra 
of isolated, strongly magnetized, thermally emitting neutron stars. 
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1. Introduction For interpretation of the XDINS spectra, it may be neces- 
sary to take the phenomenon of "magnetic condensation" into 
Recent observations of neutron stars have provided a wealth account. The strong magnetic field squeezes the electron clouds 
?-h \ of valuable information, but they have also raised many new around the nuclei a nd thereby increase s the binding and cohe- 
. questions. Particularly intriguing is the class of radio-quiet neu- sive energies (e.g.. lMedin & LailllOOrjl and references therein), 
tron stars with thermal-like spectra, commonly known as X-ray Therefore XDINSs may be "naked," with no appreciable at- 
dim isolated neutron s tars (XDINSs ), or the Magnifice nt Seven mosphere above a condensed surface, as first conjectured by 
(see, e.g., reviews by lHaberll I2007I and ITurolla 2009, and ref- IZane et al J (I2002I) . or thev may have a relatively thin atmosphere, 
erences therein). Some of them (e.g., RX J1856.5-3754) have with the spectrum of outgoing radiation affected by the proper- 
featureless spectra, whereas others (e.g., RX J1308.6+2127 and ties of th e condensed surfac e beneath the atmosphere, as sug- 
RX J0720.4-3 125) have broad absorption features with energies gested by Motc h et al.l (120031) . 

~ 0.2-2 keV. In recent years, an accumulation of observational Reflectivities of condensed metallic surfaces in strong 

evidence has suggested that XDINSs may have mag netic fields magnetic fields h ave been s t udied in several papers 

B ~ 1 0° - 10 14 G and be related to magnetars (e.g.,|Mereghettj| (iBrinkmannl 119801; ITurolla et al l 120041; Ivan Adelsberg et all 



|2008|). 200 51: iPerez-Azorm et al.l I2005D . IBrinkmann (1980) and 



Turol la et al.1 (12004ft neg lected the motion o f ions in the con- 

densed mat ter, whereas | van Ade lsberg et al. (2005) (hereafter 

e-mail: palex@ast.ro . ioffe . ru Paper I) and lPerez-Azorfn et al.l d2005h considered two opposite 
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limiting cases: one that neglects the ion motion ("fixed ions") 
and another where the ion response to the electromagnetic wave 
is treated neglecting the Coulomb interactions between the 
ions ("free ions"). A large difference between these two limits 
occurs at photon frequencies below the ion cyclotron frequency, 
but the two models lead to almost the same results at higher 
photon energies. We expect that in reality the surface spectrum 
lies between these two li mits (see Paper I for discu ssion). The 
results of Paper I and of iPerez-Azorfn etaL (2005) are similar 
to each ot her, but significantly differ from the earlier results. In 
particular, iTurolla et alj (120041) found that collisional damping 
in the condensed matter leads to a sharp cutoff in the emission 
at low photon energies, b ut such a cutoff is absent in Paper I and 
IPerez-Azorfn et alTd2005h . It is most likely that this difference 
arises from the "one- mode" descr i ption for the transmitted 
radiation adopted by ITurolla et al.l d2004l) (see Paper I for 
details). All the previous works relied on a complicated method 
of finding t he tra nsmitted radiation modes, originally due to 
Brinkmann (1980). We replace it with a more reliable method 
described below. 

IHo et alJ d2007i) (see also lHoll2007h fitted multiwavelength 
observations of RX J1856. 5-3754 with a model of a thin, mag- 
netic, partially ionized hydrogen atmosphere on top of a con- 
densed iron surface; they also disc ussed possible mechanis ms of 
creation of such a thin atmosphere. Suleim anov et alJ ( 120091) cal- 
culated various models of fully and partially ionized finite atmo- 
spheres above a condensed surface including the case of "sand- 
wich" atmospheres, composed of hydrogen and helium layers 
above a condensed surface. 

The large variety of theoretical possibilities complicates the 
modeling an d interpretation of the sp ectra. In order to facili- 
tate this task. ISuleimanov et al.l ((2010) (hereafter Paper II) sug- 
gested an approximate treatment, in which the local spectra to- 
gether with temperature and magnetic field distributions are fit- 
ted by simple analytic functions. Being flexible and fast, this 
approach is suited to constrain stellar parameters prior to per- 
forming more accurate but computationally expensive calcula- 
tions of model spectra. The reflectivity of the condensed sur- 
face was modeled by a simple steplike function, which roughly 
described the polarization-averaged reflectivity of a magnetized 
iron surface at B — 10 13 G, but depended neither on the mag- 
netic field strength B nor on the angle <p between the plane of 
incidence and the plane made by the normal to the surface and 
the magnetic field lines. 

In the present work, the numerical method of Paper I and 
the approximate treatment of Paper II are refined. We develop 
a less complicated and more stable method of calculations and 
construct more accurate fitting formulae for the reflectivities of a 
condensed, strongly magnetized iron surface, taking the depen- 
dence on arguments B and <p into account. The new fit reproduces 
the feature near the electron plasma energy, obtained numeri- 
cally in Paper I but neglected in Paper II. Two versions of the fit 
are presented in Sect.|2]for the models of free and fixed ions dis- 
cussed in Paper I. In addition to the fit for the average reflectivity, 
we present analytic approximations for each of the two polariza- 
tion modes, which allow us to calculate the polarization of radi- 
ation of a naked neutron star. In Sect. [3] we consider the radiative 
transfer problem in a finite atmosphere above the condensed sur- 
face, including the reflection from the inner atmosphere bound- 
ary with normal-mode transformations, neglected in the previous 
studies of thin atmospheres. Conclusions are given in Sect. [4] In 
Appendix [A] we describe the method of calculation for the re- 
flectivity coefficients, which is improved with respect to Paper I. 



In Appendix [B] we describe an analytic model of normal-mode 
reflectivities at the inner boundary of a thin atmosphere. 

2. Spectral properties of a strongly magnetized 
neutron star surface 

2.1. Condensed magnetized surface 

Most of the known neutron stars have magnetic fields B much 
larger than the natural atomic unit for the field strength Bq = 
e 3 mlc/h 3 = 2.35 x 10 9 G that is set by equating the electron 
cyclotron energy 



neB/m e c = 1 15.77 B 13 keV 



(1) 



to the Hartree unit of energy m e e 4 /h 2 . Here, m e is the elec- 
tron mass, e is the elementary charge, c is the speed of light 
in vacuum, fr is the Planck constant divided by 27r, and B^ = 
B/10 13 G. Fields with B » Z?o profoundly affect the properties 
of atoms, molecules , and p lasma (see, e.g. jHaensel et al.l l2007, 
chap. 4). lRuderma"nl d!971l) suggested that the strong magnetic 
field may stabilize linear molecular chains (polymers) aligned 
with the magnetic field and eventually turn the surface of a neu- 
tron star into the metallic solid state. Later studies have pro- 
vided support for this conjecture, although the surface density 
p s , and especially the critical temperature r cl -i t below which such 
condensation occurs, remain uncertain. Order-of-magnitude es- 
timates suggest 



Ps 



8.9 x 10 3 tjAZ' - 6 B\f g cm" 3 , 



(2) 



where A and Z are the atomic mass and charge numbers, and 
rj ~ 1 is an unknown nu merical fa ctor, which absorbs the 
theoretical uncertainty (see lLail 120 01). The value 77 = 1 cor- 
respon ds to the equa tion of state provided by the ion-sphere 
model dSalpeterll954l) . Results of the zero-temper ature Thomas- 
Ferm i model for S6 Fe at 10 10 G < B < 10 13 G dFushiki et al.l 
119891; iRog nvaldsso n et alJ 119931) can be approximated (within 
4%) by Eq. © with 77 * 0.2 + 0.0028/7^|_ 6 , whereas the finite- 



temperature Thomas-Fermi model of Thorolfss on et al.l ([l998) 
does not predict magnetic condensation at all. The most compre- 
hensive study of cohesive pro perties of the magnetic con densed 
surface has been conducted bv lMedin & Lai! (2006, 2007), based 
on density-functional theory (DFT). iMedin & Lail {2006) cal- 
culated cohesive energies Q s of the molecular chains and con- 
densed phases of H, He, C, and Fe in strong magnetic fields. 
A comparison with previous DFT calculations by other authors 
suggests that Q s may vary within a factor of two a t B > 10 12 G, 
depen ding on the approximations em ployed (see Medin & Lai 
2006 for references and discussion). Medi n & Lai (2007) cal- 
culated equilibrium densities of saturated vapors of He, C, and 
Fe atoms and chains above the condensed surfaces and obtained 
r cr i t at several values of B by equating the vapor de nsity to p s . 
Unlike previous authors. [Medin & Lail (2006. 2007) have taken 
the electronic band structure of the metallic phase into account 
self-consistently. However, in the gaseous phase, they still did 
not allow for atomic motion across the magnetic field and did 
not take into account a detailed tr e atmen t of excited atomic 
and molecular states. Medi n & Lail ([2007) calculated the sur- 
face density assuming that the linear molecular chains (directed 
along E) form a rectangular array in the perpendicular plane and 
that the distance between the nuclei along the field lines is the 
same in the condensed matter as in t he separate molecula r chain 
(lMedin&Laill2006t lMedinll2012t) . IMedin & Lail (120071) found 
that the critical temperature is r cr i t * 0.08(2s/^b- Their numer- 
ical results for 56 Fe at 0.5 < B13 < 100 can be described by 



2 



A. Y. Potekhin et al.: Spectra of neutron stars with metallic surfaces 



expression r cl ; t « (5 + 2 Z?n) x 10 5 K for the critical tempera- 
ture and by Eq. (|2]i with 77 w 0.55 for the surface density, with 
uncertainties below 20% for both quantities. An observational 
determination of the phase state of a neutron star surface would 
be helpful for improving the theory of matter in strong magnetic 
fields. 

The density of saturated vapor abo ve the condensed sur - 
face rapidly decreases with decreasing T dLai & Salpeter|[T997l) . 
Therefore, although the surface is hidden by an optically thick 
atmosphere at T « T^t. the atmosph ere becomes optical ly thin 
at T <K r cl j t . Also, as suggested by iMotch et al.l d2003l) . there 
may be a finite amount of light chemical elements (e.g., H) on 
top of the condensed surface of a heavier element (e.g., Fe). In 
addition, the same atmosphere may be optically thick for low 
photon energies and transparent at high energies. The energy at 
which the total optical thickness of a finite atmosphere equals 
unity depends on the atmosphere column density, which can in 
turn depend on temperature. At a fixed energy, the optical thick- 
ness of the finite atmosphere is different for different photon po- 
larizations, therefore the atmosphere can be thick for one po- 
larization mode and thin for another. One should take all these 
possibilities into account while interpreting observed spectra of 
neutron stars. 



2.2. Formation of the spectrum 

2.2.1. Normal modes and polarization vectors 



It is well known (e.g., lGinzburgl [l970) that under typical con- 
ditions (e.g., far from the resonances) electromagnetic radiation 
propagates in a magnetized plasma in the form of extraordinary 
(X) and ordinary (O) normal modes. These modes have different 
polarization vectors ex and eo, absorption and scattering coef- 
ficients, and refra c tion a nd reflection coefficients at the surface. 
iGnedin & Pavlovl (Il973l) studied conditions for the applicabil- 
ity of the normal-mode description and formulated the radiative 
transfer problem in terms of these modes. 

Fo l lowing the work s of [Shafranoy] (1 19671) and iGinz burg 
(fl^70T) . lHo&Lail(l200ll) derived convenient expressions for the 
normal mode polarization vectors in a fully ionized plasma for 
photon energies E much higher than the electron plasma energy 



£ pe = (47r/z 2 £> 2 n e /m e ) 1/2 * 0.0288 y/pZ/A keV, 



(3) 



where p is the density in g cm -3 . In the complex representation 
of plane waves with E oc ee l( *' r ~ a " ) , in the coordinate system 
where the z-axis is along the wave vector k, and the magnetic 
field B lies in the (xz) plane, the polarization vectors are 



e M (a) = 



1 

VI + \K M (a)\ 2 + \K z , M (a)\ 2 



iK M (a) \ 
1 

\K zM (a) 



(4) 



We use the notation M — X and M — O for the extraordinary and 
ordinary polarization modes, respectively. Km(o) and K z m(o) 
are functions of the angle a between B and k. They are deter- 
mined by the dielectric tensor of the plasma and thus depend 
on the phot on energy E, as w ell as p, B, T and the chemical 
composition. iHo & Laild2003i) calculated Km and studied the po- 
larization of normal modes includi ng the effect of t he ele ctron- 
positron vacuum polarization, while Potekh in et alj ([2004 ) addi- 
tionally took into account incomplete ionization of the plasma. 




Fig. 1. Illustration of notations. The z axis is chosen perpendic- 
ular to the surface, and the (xz) plane is chosen parallel to the 
magnetic field lines, which make an angle 6b with the normal to 
the surface. The direction of the reflected beam with wave vector 
kf is determined by the polar angle 6k and the azimuthal angle <p, 
and a T is the angle between the reflected beam and the field lines. 
Thick solid lines show the reflected beam and magnetic field di- 
rections, thin solid lines illustrate the coordinates, and dashed 
lines show the incident photon wave vector k t and its quadrant. 
The lines marked e^ illustrate the basic polarizations adopted 
for the description of reflectivities: el \ and el 2 are perpendicu- 



lar to the wave vectors k[ and k\, respectively; e^ are parallel 

to the surface, and e£ lie in the perpendicular plane. The axes 

x? and y lie in the plane made by ef and e\ , , x' being aligned in 
the plane made by B and k r . 



2.2.2. Emission and reflection by a condensed surface 

The condition E > E pe is usually satisfied for X-rays in neutron 
star atmospheres, but not in the condensed matter. Let us con- 
sider a surface element that is sufficiently small for the variation 
of the magnetic field strength and inclination to be neglected. 
We shall treat this small patch as plane, neglecting its curvature 
and roughness. We choose the Cartesian z axis perpendicular to 
this plane and the x axis parallel to the projection of magnetic 
field lines onto the xy plane. We denote the angle between the 
field and the z axis as 6b, the incidence angle of the radiation as 
6k, and the angle in the xy plane made by the projection of the 
wave vector as (p (Fig. [T]i. The angle between the wave vector 
and magnetic field lines is given by 



cos aj, r = sin 6b sin 6k cos <p + cos 63 cos Ok 



(5) 



for the incident and reflected waves, respectively. The surface 
emits radiation with monochromatic intensities 



I EJ =JjB E /2 (7 = 1,2). 



(6) 



Here, the basis for polarization is chosen such that the waves 
with 7 = 1 and 2 are linearly polarized parallel and perpendicular 
to the incident plane, respectively (Fig. [TJ; Igj dZ dD. dE dt gives 
the energy radiated in the jth wave by the surface element dZ, in 
the energy band (E, E + dE), during time dt, in the solid angle 
element dQ around the direction of the wave vector k. We use 
the function 



Br 



By_ 

2nti 



47r 3 H 3 c 2 (e EfkBT - 1)' 



(7) 



where B Y is Planck's spectral radiance and is the Boltzmann 
constant. Dimensionless emissivities for the two polarizations, 
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normalized to blackbody values, are Jj = 1 — Rj, where Rj is the 
effective reflectivity of mode j defined in Appendix [A] The re- 
flectivities depend on surface material, photon energy, magnetic 
field strength B and inclination 9b, and on the direction of k. 
For nonpolarized radiation, it is sufficient to consider the mean 
reflectivity R = (Ri + Ri)j2 (see Paper I). 



2.2.3. Reflectivity calculation 

A method for calculating reflectivity coefficients was devel- 
oped in Paper I. However, it is not easy to implement. Though 
it mostly produces correct results, in some ranges of model 
parameters it can yield unphysical results, which are difficult 
to distinguish from the correct ones. In the present work we 
present an improved method that avoids this complication (see 
AppendixlAl. 

Using our new method, we have calculated the spectral prop- 
erties of a condensed Fe surface and compared the results with 
those in Paper I. As in Paper I, we consider two alternative mod- 
els for the response of ions to electromagnetic waves in the con- 
densed phase: one neglects the Coulomb interactions between 
ions, while the other treats ions as frozen at their equilibrium 
positions in the Coulomb crystal (i.e., neglecting their response 
to the electromagnetic wave). 

In the first limiting case (thick lines in Fig.[2j, the reflectivity 
exhibits different behavior in three characteristic energy ranges: 
E < E c i, £ c i < E < Eq, and E > Eq, where 



E ci = hZeB/Am u c = 0.0635 (Z/A) fi B keV 



(8) 



is the ion cyclotron energy, m u is the unified atomic mass unit, 
and 



Ec - E c 



(9) 



In addition, there is suppression of the reflectivity at E ~ E l 
the exact position, width, and depth of the suppression depend 
on the geometry defined by the angles 9b, 9k, and (p. 

In the opposite case of immobile ions (thin curves in Fig. [2]), 
the reflectivity has a similar behavior at E > E c \, but differs at 
E < E c c, it does not exhibit the sharp change at E as E cl , but 
smoothly continues to the lower energies. As argued in Paper I, 
we expect that the actual reflectivity lies between these two ex- 
tremes. 

The new results, shown in Fig. [2] display the same qualita- 
tive behavior as in Paper I, but exhibit considerable deviation 
from the previous calculations for some geometric settings in 
the energy range £ c ; < E < Ec- Thus, the qualitative results and 
conclusions of Paper I are correct, but the new method described 
in Appendix [A] is quantitatively more reliable. 

If 9b = 9k — 0, then an approximate analytic solution (ne- 
glecting the finite electron relaxation rate in the medium; see 
Paper I for discussion) is R » (R + + R-)/2, where 



R 



1 



n? + 1 



„(0) 



1 ± 



Ece(E ± E c i) 



1/2 



(10) 



Compared to the numerical results, Eq. (TTOb provides a good ap- 
proximation at E < E cl . Therefore, we use it in the analytic fit 
described below. 



2.3. Results for iron surface 

In the numerical examples presented below we assume a con- 
densed 56 Fe surface and use the estimate of the surface density 
given by the ion-sphere model - that is, we set rj — 1 in Eq. (0. 
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Fig. 2. Dimensionless emissivity J = 1 - R as a function of pho- 
ton energy E for a condensed Fe surface with B = 10 13 G and 
T = 10 6 K. The top panel shows several cases with varying 9b 
and fixed 9k = n/4, <p = n/4. The bottom panel shows several 
cases with varying <p and fixed 9k - n/4, 9b - n/4. These plots 
should be compared with Figures 5 and 6 of Paper I. 

2.3.1. Mean reflectivity 

In practice, the average normalized emissivity / = 1 - R is usu- 
ally more important than the specific emissivities Rj. In Paper II, 
R(E) was replaced by a constant in each of the three ranges men- 
tioned in Sect. 12.2. 3l (£' < E c u E c \ < E < Eq ~ Eq, and E > Ec). 
For simplicity, the values of these three constants were assumed 
to depend only on 9b and 9k, but not on <p or B. Here we propose 
a more elaborate and accurate fit, which is a function of E, B, 9b, 
9k, and tp, for the magnetic field range 10 12 G < B < 10 14 G and 
photon energy range 1 eV < E < 10 keV. In the approximation 
of free ions, the average reflectivity of the metallic iron surface 
is approximately reproduced by 



J = 



J & in Region I, 



Jb(1-Jc) + 



Jc_ 
1 +L 



in Region II. 



(11) 



Region I is the low-energy region defined by the conditions 
E < E c [ and 7a > Jb- Region II is the supplemental range of 
relatively high energies in which either of these conditions is vi- 
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Fig. 3. Dimensionless emissivity / = 1 - R as a function of 
photon energy E for condensed Fe surface at B — 10 13 G (top 
panel), B = 10 12 G and 10 14 G (bottom panel), with magnetic 
field lines normal to the surface, for two angles of incidence ft = 
and ft = tt/4, as marked near the curves. Solid lines show 
our numerical results, and dashed lines demonstrate the fit. For 
comparison, dotted lines reproduce the simplified approximation 
used in Paper EL 

olated. The functions J a, Jb, and Jc are mainly responsible for 
the behavior of the emissivity at E < E C i, E c \ < E < Ec, and 
E > Ec, respectively, while the function L describes the line at 
E * £p e . The value 



E c = E a + £^ e /£ ce (12) 
is the energy at which the square of the effective refraction index 



h 2 = \- 



E ce (E E cl ) 

(analogous to Eq. ( fTUl i) becomes positive with increasing E in 
the range E > E cl . In Eqs. ( fT2l and ( fT3l >. 



y3 - 2 cos ft. 
The low-energy part of the fit in Eq. (fTTT i is given by 
J A = [l-A(E)]J (E), 
where 



(14) 



(15) 



= ^-cosft + 
2 VI +B l3 



0.7- 



0.45 
/o(0) 



(sin ft) 4 (1 -cos a), (16) 



Jo(E) = 1 - ±(fl™ + Rf> ), and Rf> are given by Eq. (T0]>. 

Accordingly, / (0) = 4 ( ^JEdE^i + l)" 1 ( ^jE~IE~c~ + l)"'. In 
Eq. ( fTol l and hereafter, <? without subscripts denotes min(a r , af). 



In the intermediate energy range, E cl < E < Ec, there is a 
wide suppression of the emissivity. We describe this part by the 
power-law interpolation between the values at E c i and Ec' 



J B = (E/E c yj(E c ), where p = 



ln[J(E c )/J(E ci )] 
ln(£ c /£ci) 



(17) 



The values J(E C1 ) and J(Ec) are approximated as follows: 



J(Ec) 



J(E ci ) = 



1 



0.05 



-(1 + cos 8b sin ft) 



2 l+B 13 
-0.15(1 - cos ft) sin a, 



2«o 



(1+rco) 2 



1 + 



cos 6b - cos ft 
2(1 +B a ) 



(18) 
(19) 



where n = (1 + £ 2 e /2£ ce £ ci ) I/2 . 

The steep slope at £ > £c is described by Eq. ( fTOb with fipe 
replaced by E pe : 

|4n/(l+n) 2 at£>£ c , 
\ at E < £ c , 

« being given by Eq. ( fP3l ). 

Finally, the lowering of /(£) at E > Ec is fit by 



(20) 



Q.lTEpJEc 



1 +X 4 



+ 0.21e- (£/£ >" )2 



(sinft) 2 W L 



(21) 



^ g 

x = L (i -cosftr 1 



2£ pe W L 

£ L = £ pe [l + 1.2(1 - cos ft ) 3/2 ] [l - (sin 6 B ) 2 /3] , 

W L = 0.8 (E c /E pe f 2 Vsin(ff/2) [l + (sin 0b) 2 ] . 

The line at Ei disappears from the fit (L — > 0) when radiation 
is parallel to the magnetic field (a — > 0). This property is not 
exact; our numerical results reveal a remnant of the line at a — > 
0, which is relatively weak, but may become appreciable if the 
magnetic field is strongly inclined (dg > tt/4). 

Examples of the numerical results for the normalized emis- 
sivities are compared to the analytic approximation in Figs. [3]- 
[5] For most geometric settings, the fit error lies within 10% in 
more than 95% of the interval -3 < log I0 £'(keV) < 1. The 
remaining < 5% are the narrow ranges of E where J(E) sharply 
changes. Exceptions occur for strongly inclined fields (6b > tt/4) 
and small <p, where the error may exceed 20% in up to 10% of 
the logarithmic energy range. 

Our fit does not take into account the dependence of re- 
flectivity on temperature. Temperature of the condensed matter 
enters in the calculations through the effective relaxation time, 
which determines the damping factor (see Paper I). This dis- 
regard is justified by the weakness of the T-dependence of the 
results. With decreasing T, the transitions of R(E) between char- 
acteristic energy ranges become sharper, and the feature at Ei 
becomes stronger. The bulk of our calculations employed in the 
fitting was done at T ~ 10 6 K. In Fig. [4] where the solid lines 
represent the numerical results at T = 10 6 K, an additional dot- 
ted line is drawn at T = 3 x 10 5 K, in order to illustrate the 
T-dependence. 

A small modification of the proposed approximation can de- 
scribe the alternative model of fixed ions (see Paper I). In this 
case, it is sufficient to formally set E Q i — > in the above equa- 
tions and replace Eq. (TTTb by 



Jb = 



J(Ec) 



l-p+p(E c /Ey 



0.6' 



where p = 0.1 



1 + sin ( 
1 +B, 



(22) 
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Fig. 4. Dimensionless emissivity J — 1 - R as a function of 
photon energy E for a condensed Fe surface with inclined mag- 
netic field (B — 10 13 G, 6b = and inclined incidence of 
radiation (0* = n/4 and four values of <p listed in the figure). 
Solid lines show our numerical results for the model of free ions 
at T = 10 6 K, and short-dashed lines demonstrate the fit. For 
comparison, the dotted line reproduces our numerical results for 
(f = n/2 and T = 3 x 10 5 K. 




0.001 



Fig. 5. Same as in Fig. |4] but for B = 10 14 G and two values of 
<p. For comparison, the dotted and long-dashed lines reproduce 
our numerical results and analytic approximation, respectively, 
for the model of fixed ions. 



As an example, in Fig. the dotted lines illustrate the numeri- 
cal results in the model of fixed ions, and the long-dashed lines 
correspond to Eq. (l22l >. 

2.3.2. Polarization 

We also constructed approximations for the emissivities in each 
of the two modes. Their functional dependence on E and geo- 
metric angles in Fig. Q] is more complicated than the analogous 
dependence for the average J = (Ji +J2)/2. We do not accurately 
reproduce these complications, in order to keep the fit relatively 
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Fig. 6. Dimensionless emissivity for the linear polarization e\, 
J i = 1 - Ri, as a function of photon energy E for condensed 
Fe surface at different magnetic field strengths and geometric 
settings. Top panel: magnetic field B = 10 13 G is normal to the 
surface, and angles of radiation incidence are 6^ = 0, 7r/6, tt/4, 
and 7r/3. Bottom panel: B - 10 13 5 G, magnetic field lines and 
the photon beam are both inclined at Be = 0k = tt/4, and the 
azimuthal angle takes values <p — 0, 7r/4, and n/2. Solid lines 
show the numerical results, and dashed lines demonstrate the fit. 



simple, but reproduce general trends. For j = 1, our free-ion 
approximation has the same form as Eq. (fTTT i: 



h = 



Jai in Region I, 

Jbi(1-Jc) + MI-Rl) in Region H. 



(23) 



Here, we retain 7c given by Eq. d20l >. The shape of the line near 
the plasma frequency is unchanged and is described by Eq. (|2H . 
but the line strength is different, since L enters Eq. d23l through 
the function 



R L = (sin0 B ) 1/4 [2 -(sin a) 4 ] j^. 



(24) 



The functions that describe emissivity in mode 1 at E < Eq are 
Jai = [1 -Ai]/ A , (25) 



Cl\ 



1 +0.6Bi 3 (cos6» B ) 2 ' 



fli = 1 — (cos 8b) 2 cos 6k - (sin 6b) 2 cos a; 
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Fig. 7. Degree of linear polarization P\ m [Eq. d28l ll as a function 
of photon energy E for condensed Fe surface. The values of B, 
directions of the field and the wave vector, and line types are 
same as in Fig. [6] 



7bi = (E/E C ) P1 ME C ), Px = 



\n\J\(E c )Ux(E^\ 
ln(E c /E cl ) " 



(26) 



J 1 (E ci ) = (l-a l )J(E ci ), 

~ 1 0.05 sin6 B 

JliEc) = 2 + TT^I + — 

In the fixed-ion case, it is sufficient to set E c 
Eq. (ESJ by 



Ji(Ec) 



0.1 +0.9 (Ed E)' 



0.4 



and to replace 



(27) 



For the second mode, no additional fitting is needed, because 
R 2 = 2R - Ri and J 2 = 2J - J x . 

Figure [6] compares the use of Eqs. (l23l - d26li (dashed lines) 
to numerical results (solid lines). The upper panel shows the case 
where the field lines are perpendicular to the surface. In this case 
the line at Ei disappears from mode 1 . Hence the line in R seen 
in Fig.[3]for k + is entirely due to the mode 2. As soon as the 
field is inclined, the line is redistributed between the two modes 
(the lower panel of Fig. [(J. In the latter case the numerical results 
(solid lines) show a more complex functional dependence R\(E) 
in the range E c \ < E < Eq, which is not fully reproduced by our 
fit (dashed lines), for the reasons discussed above. 

Note that the azimuthal angle tp enters the fit only through a. 
As a consequence, the fit is symmetric with respect to a change 
of sign of tp. This property may seem natural at the first glance, 
however, we note that the numerical results do not strictly obey 



this symmetry, which holds for the nonpolarized beam, but not 
for each of the polarization modes separately. We have checked 
that this is not a numerical artifact: because the magnetic field 
vector B is axial, there is no strict symmetry with respect to 
the (x, z) plane. A reflection about this plane would require si- 
multaneous inversion of the B direction in order to restore the 
original results. However, as long as the electromagnetic waves 
are nearly transverse (i.e., K z in <j4j and ( IA.24b are small), the 
asymmetry is weak, allowing us to ignore it and thus keep the fit 
relatively simple. 

The analytic approximations in Eqs. ( fTTT i and d23l allow one 
to evaluate the degree of linear polarization of the emitted radia- 
tion 



Pun = (Ji - Ji)/2J = (R 2 - Ri)/(2 - 2R). 



(28) 



For example, the two panels of Fig. [7] show Pi; n for the same 
directions of the magnetic field and the photon beam as in the 
respective panels of Fig. [6] We see that the analytic formulae, 
originally devised to reproduce the normalized emissivities, also 
reproduce the basic features of Pii n (£'). Although the feature at 
E ~ £p e is absent in the top panel of Fig. |6l it reappears in the 
top panel of Fig.|7]due to the contribution of R 2 in Eq. ( |28| >. 



3. X-ray spectra of thin atmospheres 

3. 1 . Inner boundary conditions 

Propagation of radiation in an atmosphere is described by two 
normal modes (see Sect. l2.2.TT i. At the inner boundary of a thin 
atmosphere, an incident X-mode beam of intensity I* gives rise 
to reflected beams in both modes, whose intensities are propor- 
tional to and analogously for an incident O-mode. Therefore, 
the inner boundary conditions for radiation transfer in an atmo- 
sphere of a finite thickness above the condensed surface can be 
written as 



hMOk, <P)Be(T) + Rxx(0 k , <p) If(n - k , if) 



+R xo (6 k ,ip)I { £(n. 



ixxt 

e k ,<p), 



(29) 
(30) 



I%(e k ,<p) = \Jo(e k ,ip)B E (T) + R oo (e k ,<p)I°(n-0 k ^) 

+Rox(0k,<p)I%(n-0k,<p), 

where 1% (M = X, O) are the specific intensities of the X- and 
O-modes in the atmosphere at p — p s , Rmm 1 are coefficients 
of reflection with allowance for transformation of the incident 
mode M' into the reflected mode M, and Jm are the normalized 
emissivities. The latter can be written by analogy with as 
Jx = 1 — Rx an d Jo — 1 - ^o, where Rx = Rxx + Rxo and 
#o = floo + Rox (c f. Paper I). 

iHo et al.l d2007l) retained only the emission terms ^JmBe on 
the right-hand sides of Eqs. ( l29l . ( l30l l. The reflection was taken 
into account in Paper II, but calculations were performed ne- 
glecting Roo, Rox, and Rxo, under the assumption that Rxx is 
equal to R and does not depend on <p. Here we use a more re- 
alistic, albeit still approximate, model for Rmm', described in 
Appendix IE1 



3.2. Results 

Here, we illustrate the importance of the correct description of 
the reflection for computations of thin model atmospheres above 
a condensed surface. To this end, we have calculated a few model 
atmospheres with normal magnetic field (therefore, 0b — <p = 0, 
and ff r = QTj = k ), taking the model with B = 4 x 10 13 G, 
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Fig. 8. Emergent spectra (top panel) and temperature struc- 
tures (bottom panel) for the fiducial model atmosphere (solid 
curve) and for model atmospheres which are calculated using the 
fixed-ions approximation for the reflectivity calculations (dashed 
curves), and the inner boundary condition from Paper II (dotted 
curves). In the top panel the diluted blackbody spectrum which 
fits the high-energy part of the fiducial model spectrum is also 
shown (dash-dotted curve). 



effective temperature - 1.2 xlO 6 K, and surface density 2 = 
10 g cirr 2 as a fiducial model. In the fiducial model the free-ions 
assumption for condensed-surface reflectivity is used. 

For these compu t ations we use the numerical code described 
in Suleimanov et al . (2009), with a modified iterative procedure 
for temperature corrections. We evaluate thes e corrections using 
the Unsold-Lucy method (e.g. jMihalasl 1978b . which gives a bet- 
ter convergence for thin-atmosphere models than other standard 
methods. In our case, the deepest atmosphere point is the up- 
per point of the condensed surface. The temperature correction 
at this point is obtained as follows: the total flux at the bound- 
ary between the atmosphere and the condensed surface is fixed 
and, therefore, the following energy balance condition has to be 
satisfied: 



4tt 2 
= B tot k Rh + JR + H-. 



i r°° r 1 

-J dE j i (/IOu) + pd, 



(31) 



Here, cr SB is the Stefan-Boltzmann constant, ji = cos 6^, and 

Jr»co 
B E dE , 
o 
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Fig. 9. Top panel: dimensionless emissivities for coefficients of 
reflection Rxx (dashed curve), Rxo (dotted curve), Rox (dash- 
dot-dotted curve), and Roo (dash-dotted curve). The quantities 
are calculated at the bottom of the fiducial model atmosphere 
for angle between the radiation propagation and magnetic field, 
9k = 10°, together with the total dimensionless emissivity (solid 
curve). Bottom panel: dimensionless outward specific intensities 
(inner boundary condition) at the bottom of the fiducial model 
atmosphere for the X-mode (solid curve) and O-mode (dashed 
curve). For comparison, the dotted curve shows the same for 
the X-mode, calculated using the inner boundary condition from 
Paper II (in this case the dimensionless specific intensity of the 
O-mode equals 0.5). 



km, = 



2B U 



JR 



-|— f B E dE f (l-R) f id f i, 

K tot Jo Jo 

i r°° 

- 2 J. " £ 
x J (/fOu)(«xx+/?ox) + /°(M)(«xo+«oo))yudyu ( 

-I pCO M) 

H ~ = 2j dE J i ( / ^) + /°(/"))^d / z. (32) 

Generally, the condition OTb is not fulfilled at a given tem- 
perature iteration. Therefore, we perform a linear expansion of 
the integrated blackbody intensity: 



Ho = (B,ot + AB tot )£ RL + JR + H, 

and find a corresponding temperature correction 



(33) 



(34) 
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This last-point correction procedure is stable and has a conver- 
gence rate similar to the Unsold-Lucy procedure at other depths. 

We have also changed the depth grid for a better descrip- 
tion of the temperature structure in thin-atmosphere models. In 
semi-infinite model atmospheres that do not have a condensed 
surface as a lower boundary, a logarithmically equidistant set of 
depths is used. However, in thin-atmosphere models, such a set 
yields insufficient accuracy at the boundary between the atmo- 
sphere and condensed surface. To improve the description of the 
boundary, we divide the model atmosphere into two parts with 
equal thicknesses and use logarithmically equidistant depth grids 
for each of them. In the upper part the grid starts from outside 
(the closest points are at the smallest depths), while in the lower 
part it starts from the condensed surface (the closest points are 
at the deepest depths). This combined grid allows us to describe 
the whole atmosphere with the desired accuracy of 1 % for the 
integral flux conservation. 

In Fig. [8] we show the emergent spectra and temperature 
structures for three different model atmospheres with the same 
fiducial set of physical parameters. The model spectrum com- 
puted using the inner boundary condition described in Paper II 
(the "old model") significantly differs from the two other model 
spectra computed using the improved boundary conditions of 
Eqs. d29j >, (f3Qb . The latter two models are calculated using the 
free (fiducial model) and fixed ions (alternative model) assump- 
tions for the condensed surface reflectivities. 

The differences between the spectra and the temperature 
structures of these two models are very small. The atmo- 
sphere temperature near the condensed surface with fixed ions is 
slightly smaller than the temperature near the condensed surface 




E (keV) 

Fig. 11. Top panel: total emergent spectrum of the fiducial 
model (solid curve) together with emergent spectra in the X- 
mode (dashed curve) and in the O-mode (dash-dotted curve). 
The blackbody spectrum with T = is also shown (dotted 
curve). Bottom panel: emergent specific intensities of the fidu- 
cial model for six angles 9. 



with free ions. The flux in the spectrum of the fiducial model is 
approximately twice that of the alternative model at photon en- 
ergies E smaller than the iron cyclotron energy E c \ = 0.118 keV. 
At larger energies the spectra are very close to each other. We 
note that the old model has been computed using the free-ions 
assumption. 

The difference in the emergent spectra between the old and 
new model atmospheres is significant. In the old model, there is 
a deep depression of the spectrum between E c i and Ec with an 
emission-like feature around the absorption line at the proton cy- 
clotron energy E cp = 0.252 keV. In the new spectra this complex 
feature between Ed and Eq is completely different. The total de- 
pression is not significant, but instead of the flux increase, there 
appears a deep absorption feature at photon E > E cp . This ab- 
sorption corresponds to the bound-bound transitions in hydrogen 
atoms in strong magnetic field (Hb-b feature). 

It is clear that this difference arises due to the inner bound- 
ary condition. The bottom panels of Figs.l9land[T0lillustrate the 
difference in the outgoing flux at the boundary between the at- 
mosphere and the condensed surface for the old and new models. 
This difference is especially large for the flux in the X-mode. In 
the old model we assumed complete reflection in the X-mode; 
therefore, the reflected flux was small as the atmosphere was 
optically thin at these energies. As a result we found a small 
emergent flux at these energies. In the new models, we have sig- 
nificant mode transformation due to reflection, which causes an 
appreciable part of the energy from the O-mode to convert into 
the X-mode; the converted photons then almost freely escape 
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Fig. 13. Comparison of emergent spectra (top panel) and tem- 
perature structures (bottom panel) of the fiducial model (solid 
curves) and the model atmosphere with different magnetic field 
(B = 10 14 G), but with the same surface density (dashed curves). 



from the atmosphere. The reflectivity coefficients Rmw for two 
angles 6k are shown in the upper panels of Figs.l9land[T0l 

The total equivalent widths (EWs) of the complex absorption 
features in the spectra of the new models are smaller than the EW 
of this feature in the spectrum of the old model. Nevertheless, 
they are still significant, with EW w 220-250 eV, if the con- 
tinuum is assumed to be a diluted blackbody spectrum that fits 
the high-energy tail of the model. The parameters of the diluted 
blackbody spectrum are a color correction factor f c - T/T e g = 
1.2, and a dilution factor D = l.l -4 . The range of values for 
EW is sufficient to expla in the observed absorption features of 
XDINSs (for reviews, see lHaberlll2007l:lTurollal2609h . 

The new and old spectra are strongly polarized, with most of 
the energy radiated in the X-mode (see Fig.[TT] upper panel). We 
note that, for the parameters of the fiducial model, the vacuum 
resonance density occurs between the X and O mode photon de- 
coupling densities. Therefore, the polarization signal does not 
exhibit a rotation of the plane of polarization between low and 
high en ergies. In contras t, mo dels that exhibit this e f fect, c onsid- 
ered bv lLai & H o (2003) and lvan Adelsberg & LaH (120061) . have 
a lower magnetic field and higher effective temperature, caus- 
ing the vacuum resonance to occur outside the X and O photo- 
spheres. 



The angular distribution of the emergent flux is different in 
the two models (Fig.QTl bottom panel), especially at photon en- 
ergies between Eq and AEq- In the old model, the angular distri- 
bution is peaked around the surface normal. In the new model, 
the emergent radiation is almost isotropic, with a peak around 
the surface normal at the broad Hb-b absorption feature. 

The influence of the atmosphere thickness on its emergent 
spectrum and the temperature structure is illustrated in Fig. [12] 
A thinner atmosphere with I = I g cm" 2 has an insignificant 
Hb-b absorption feature because it is formed at higher column 
densities (» 1 - 2 g cm" 2 ). The fiducial model has the smallest 
temperature at these column densities among all the models. As 
a result, the Hb-b absorption feature is most significant in the 
spectrum of this model. The spectrum of the semi-infinite atmo- 
sphere has a hard tail and does not have any feature at E c \. 

The importance of the Hb-b absorption feature decreases if it 
is located far from the maximum of the spectrum. This is illus- 
trated by the comparison of the fiducial model with the model 
calculated for B = 10 14 G (Fig.Qj). In the latter case, the Hb-b 
absorption feature is less visible and cools the atmosphere at col- 
umn densities about a few g cm" 2 less efficiently, although the 
EW decreases insignificantly. 
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3.3. Discussion: toward models of observed spectra 

Our calculations are presented for a local patch of the neutron- 
star surface with particular values of T e ff and B. By taking sur- 
face distributions of and B into account, one can construct 
an emission spectrum from the entire neutron star; however this 
spectrum is necessarily model-dependent, as the and B dis- 
tributions are generally unknown. If these distributions are suffi- 
ciently smooth, then integration over the surface makes absorp- 
tion features broader and shallower, as demonstrated, e.g., in the 
case of cooli ng neutron star s with dipole magnetic fields and 
semi-infinite dHo et al.l l2008) or thin (Paper II) partially ionized 
hydrogen atmospheres. As shown in Paper II, smearing of the 
features is stronger, if the crustal magnetic field has a strong 
toroidal component, but weaker, if radiation is formed at small 
hotspots on the surface, whe re B can be considered a s constant. 
Using the results of Paper II, Ha mbarvan et alJ (1201 ll) fitted ob- 
served phase-resolved spectra of XDINS RBS 1223 and derived 
constraints on temperature and magnetic field strength and distri- 
bution in the X-ray emitting areas, their geometry, and the grav- 
itational redshift at the surface. The present, more detailed ap- 
proximations for the reflectivities can be directly used to refine 
these fits and constraints. 

In our numerical examples presented above, we evaluated 
the density of the condensed matter using Eq. (|2]i with 7=1. 
An eventual correction to this approximation is rather straight- 
forward, once p s is accurately known. Indeed, the density enters 
calculations through E pe cc y]rj [Eq. (0] and through the damp- 
ing factor (Paper I). The latter dependence is relatively weak, 
therefore, it is sufficient to correct E pe in the expressions pre- 
sented in Sect. [2] 

As mentioned in Sect. 13.21 the model spectra are highly po- 
larized at the stellar surface. However, the observed polariza- 
tion signal is affected by propagation of the photons through 
the neutron-star magnetosphere. In addition to redshift and light 
bending effects near the stellar surface, the mode eigenvectors 
evolve adiabatically along with the direction of the changing 
magn etic field in the magnetosphere (see, e.g., iHevl & Shavivl 
2002). The adiabatic evolution continues until the photons near 
the polarization limiting radius, r p i, which is typically many stel- 
lar radii from the surface. At r p i, the modes couple, with the 
intensities and eigenvectors frozen thereafter (in addition, sig- 
nificant circular polarization can be generated in some cases, 
for example, in radiation fr om rapidly rotating neutron stars; see 
van Adelsberg & Lai 2006, and the references therein). 

The main effect of adiabatic photon propagation in the mag- 
netosphere is on the synthetic polarization signal from a finite 
region of the neutron-star surface. (There is an additional ef- 
fect for photons propagating through a quasi-tangential region 
of magnetic field near the stellar surface; i n the majority of ca ses 
this "QT effect" can be neglected - see IWang & Laill2009l for 
details). Since the mode properties are fixed at large distances 
from the star, where the magnetic field is aligned for photon tra- 
jectories from different areas on the star, the polarization signal 
is not as diminished due to variation in the surface magnetic field 
as mi ght be expected if vacuum polarizatio n effects are ignored 
(see lHevl & Shavivll2002t IHevl et al]l2003t) . Thus, it is possible 
that polarization features of the local thin atmosphere models de- 
scribed above may be retained in spectra from a finite region of 
the neutron star. Observed spectra and polarization signals have 
been presented in the literature, employing s everal atmosphere 
models for emission from the entire surface dHevl et aljr2 003) 
and from a finite sized hotspot (Ivan Ad elsberg & Perna 20091). 



4. Conclusions 

We have improved the method of Paper I for calculating spec- 
tral properties of condensed magnetized surfaces. Using the im- 
proved method, we have calculated a representative set of reflec- 
tivities of a metallic iron surface for the magnetic field strengths 
B - 10 12 G- 10 14 G. Based on these calculations, we have con- 
structed analytic expressions for emissivities of the magnetized 
condensed surface in the two normal modes as functions of five 
arguments: energy of the emitted X-ray photon E, field strength 
B, field inclination 9b, and the two angles that determine the 
photon direction. We have considered the alternative limiting ap- 
proximations of free and fixed ions for calculating the condensed 
surface reflectivity. 

We have improved the inner boundary conditions for the ra- 
diation transfer equation in a thin atmosphere above a condensed 
surface. The new boundary condition accounts for transforma- 
tion of normal modes into each other due to reflection from the 
condensed surface. To implement this condition we suggest a 
method for calculating reflectivities Rmiw in the normal modes 
used for model atmosphere calculations, based on analytic ap- 
proximations to the reflectivities. 

We computed a few models of thin, partially ionized hydro- 
gen atmospheres to investigate the influence of the new bound- 
ary condition on their emergent spectra and temperature struc- 
tures. The allowance for mode transformations makes the com- 
plex absorption feature between E c \ and Ec less significant and 
the atomic absorption feature more important. Nevertheless, the 
equivalent widths of this complex absorption feature in the emer- 
gent spectra are still significant (« 200-250 eV) and sufficient 
to explain the observed absorption features in the spectra of 
XDINSs. Models of thin atmospheres with inclined magnetic 
fields are necessary for detailed descriptions of their spectra. We 
plan to compute such models with vacuum polarization and par- 
tial mode conversion in a future paper. 
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Appendix A: Improved reflectivity calculation 

In this Appendix, we describe several improvements to the meth- 
ods of Paper I that have enabled us to produce a general, efficient 
code, free of numerical difficulties, that computes the correct 
value of the reflectivities over the full range of parameters used 
in neutron-star atmosphere modeling. 

In general, each incoming linearly polarized wave E{ = 
Mie, and = is partially reflected, giving rise to re- 



'1 *"2 _ "^ B 2 

fleeted and transmitted field^j 



E {I) 

j 



E {t) 

j 



(A.l) 



As shown in Paper I, the dimensionless emissivities for the two 
orthogonal linear polarizations are Jj = 1 - Rj, where 



tf/ = k/in + M 0'=i>2). 



(A.2) 



The reflected field amplitudes r\\,r\%,r%\ and r 2 2 were cal- 
culated in Paper I using an eighth-order polynomial in the re- 
fraction index n ; to determine the properties of the transmit- 
ted modes. The transmitted wave can be described by two nor- 
mal modes, thus, most of the roots obtained from that polyno- 
mial represent unphysical solutions to the equations. The condi- 
tions to identify the correct roots were derived in Appendix B of 
Paper I, with the requirements that the corresponding reflectivi- 
ties satisfy R\,Ri < 1 and that the function R(E) be continuous. 
However, for some values of the model parameters, the unphys- 
ical roots can satisfy the physical constraints on the solution. 

Here we propose an improved method, based on a fourth or- 
der polynomial, which allows for easy elimination of the un- 
physical roots. 



1 We use the inverse order of the subscripts in r m j and t„,j with respect 
to that used in Paper I. 



A. 1. Transmitted modes 

A significant simplification of the equations describing the trans- 
mitted mode properties can be obtained by writing the transmit- 
ted wave vector as 



tij = —kj = rij [sin 0j (cos <pii + sin ip y) + cos 8j zl 
= sin 6k (cos tp x + sin ip y) + n z j z, 



(A.3) 



where x, y, z are unit vectors along the x, y, z axes, respec- 
tively, and the quantities cos 8j, sin di are complex numbers sat- 
isfying the condition: cos 2 0j + sin 2 6; = 1 (cf. Appendix B of 
Paper I). The second equality in Eq. dA.3b follows from Snell's 
Law, rtj sin 6j = sin 8k, and the definition n z j = rtj cos 8j. 
From Maxwell's equations for the transmitted modes, 

A E = 0, (A.4) 
A = e + rtj ® rij - nj I, (A.5) 

where e is the dielectric tensor of the medium (see Eq. (13) of 
Paper I), I is the unit tensor, and E is the electric field vector. 
If we note that n 2 = n 2 (sin 2 6j + cos 2 6 A = sin 2 k + n 2 ., and 
apply the condition det A — to obtain a nontrivial solution to 
Eq. ( IA.4I ). the result is a fourth order polynomial in n z j: 

(A.6) 

(A.7) 

(A.8) 
(A.9) 
(A. 10) 
(A.ll) 



o-Afi^j + a-in^j + fl2«? j + a\n z j + ao = 0, 
«4 = 1 + sin 2 6k I (e zz - sin 2 9k) , 
03-2 sin &k cos <p e xz / (e zz - sin 2 &kj , 

«: = '/«.}'•... + nyyjxx ~ Vxyjyx ~ Vyxjxy ~ PxyPyx, 
a I = VyyPxx ~ VxyPyx ~ VyxPxy, 
ao = HxxVyy ~ VxyVyx, 
where the coefficients have the values: 

rjxx = ex t -sin 2 ^sin 2 (^-e 2 ,/(e, z -sin 2 t ), (A.12) 

t]xy - £xv + sin 2 6k sin tp cos <p + e xz e yz / (e zz - sin 2 6k J , (A. 1 3) 

lyx = ~£xy + sin 2 6k sin ip cos tp — e xz e yz l (e^ — sin 2 9%) , (A. 14) 

rjyy = e yy - sin 2 & k cos 2 tp + e 2 _/ [e zz - sin 2 8k) , (A. 15) 

fi xx = -2 sin 6k cos ipe xz / (e zz - sin 2 6k) , (A. 16) 

Pxy - (^yz ^^*s ip — € xz sin ip) sin 8k I \ €zz ~ sin 2 8k) , (A. 17) 

Pyx = ~ (e.i- sin ip + e yz cos tp) sin 8 k / (e zz - sin 2 8k), (A. 1 8) 

7 xx = -1 - sin 2 8 k cos 2 <p/ (e z , - sin 2 8 k ) , (A. 19) 

7xy = Jyx = _ sin 2 8k sin ip cos (p/ [e zz - sin 2 6k) , (A. 20) 

7yy = -1 - sin 2 8k sin 2 tp/ [e zz - sin 2 8 k ) , (A.21) 

The fourth order polynomial defined by equations ( 1A.61 >- 

(IA.2U has much better numerical properties than the eighth- 
order polynomial described by Eq. (A4) of Paper I. We find that 
a stable, efficient method for solving Eq. ( IA.6I > can be obtained 
by defining the matrix 



M = 



' —123/(34 


-fl 2 /l24 


-fll/fl4 


-a /a4 ' 


1 














1 








I 





1 


J 



(A.22) 



and noting that the eigenvalues of M are equal to the roots of 
Eq. ( lAToT l. We use the ZGE EV subroutine of the LAPACK li- 
brarv (lAnderson et al.lll999h to compute the eigenvalues of M. 
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Of the resulting four roots, only two correspond to physical so- 
lutions for the transmitted waves. To identify the correct roots, 
we write the spatial variation of the transmitted electric field 
as Ej (r) oc exp (iaj rij ■ r/cj. The amplitude of the electric field 
must decay in the transmitted wave region, leading to the condi- 
tion: 



SmQiJ) < 0. 



(A.23) 



At all energies and angles for the range of magnetic fields B = 
10 12 - 10 15 G, condition ( lA.23t identifies two physical solutions 
to Eq. ( IA.6b . Using these values for rtj , we write the polarization 
vectors for the transmitted wave as: 



1 



1 + \K m \ 2 + \K[ 



( K® ^ 

j 

1 

jAt) 
v K z,j t 



K' 



(t) 



7xyn 2 zJ + /3 xy n zJ + a x 



>', »»;'., + B xx n zJ + a xx ' 
e xz - sin 8/c sin <pn z j - \ e xi + sin Ok cos (pn z jj K® 



e zz - sin Ok 



(A.24) 

(A.25) 
(A.26) 



A.2. Reflectivity calculation 



Once the quantities nj and ey are known, the reflectivity of the 
medium can be calculated using the boundary conditions for 
Maxwell's equations at the condensed matter surface: 



Mxz = 0, 
ABxi = 0, 



(A.27) 
(A.28) 



where AE = E® + E (r) - E (t) and AB = + B (r) - B m are 
the differences between the fields above (incident and reflected) 
and below (transmitted) the condensed surface. For the detailed 
forms of the fields, see §3.1 of Paper I. Writing out the compo- 
nents of ( IA.27t and ( IA.281 > for the two orthogonal linear polar- 
izations of the incident wave yields a system of equations for the 
amplitudes of the reflected and transmitted modes, analogous to 
Eq. (A6) of Paper I. This set of equations can be solved as two 
independent linear systems with complex coefficients, such that: 



(A.29) 



'm 






-cosy; 


cos Ok sin <p s 








- sinyj 


- cos Ok cos <p 




tl2 




- cos Ok sin <p 


- cos <p 


y. hi 


hi , 




, COS Ok cos <p 


-sirup 



where 



C = 



cos ip 

simp 
- cos 0k sin <p 
cos Ok cos tp 



cos Ok sin <p 
- COS Ok cos tp 
COS ip 

siny> 



-1 

e (t) 

■ n (t) (0 

sin Ok cos tpe\ - n z \e\ ' 
sin Ok sin tpef^ - n z j 



-1, 

.(0 
2..r' 
« 



e 



sin Ok cos <pe 2 - n z ^e 2 \ 



sin Ok sin ipe® - n z ^ 



ande® = e®-x, e^ z 



ef-z. We solve the complex systems using 
the Z GESV subroutine of the LAPACK library (lAnderson et alJ 
Il999h . 

The corrected results for the case of a magnetized iron sur- 
face (to be compared with Paper I) are presented in Sect. 12.2.31 



Appendix B: Approximations for reflectivities at the 
bottom of a thin atmosphere 

B.1. Reflectivities of the normal modes in terms of r mj 

In general, the interface between the thin atmosphere and mag- 
netic condensed surface has reflection and transmission proper- 
ties different from those of the condensed surface in vacuum. 
Therefore, a separate calculation of the reflectivity coefficients 
r m j is needed for every set of atmosphere parameters. However, 
assuming that the atmosphere is sufficiently rarefied, we may 
approximately replace these coefficients by those in vacuum. 
Under these conditions, the plane waves in the atmosphere are 
almost transverse, hence we can approximately set K Z) m — > in 
Eq. ©. Then each incident and reflected wave can be expanded 
over the linear polarization vectors ei and ei that have been em- 
ployed in the reflectivity calculation. 

For the incident (i) and reflected (r) beams, we define or- 



thonormal vectors e 



(i,r) 
1 



z x k/\i x k\ = z x k;, r /| sinOkl 



if, and ef = ef 
vector along k. In the notations of Fig. Q] 



ev — ki x ey, and ey - ey x k r , where k denotes the unit 



ki, r = sin Ok cos <p x + sin Ok sin <pf T cos Ok z, 



(i) _ 
l 

Cv) 

2 



£j = - sin (/>x + cos ip y, 

- cos Ok (cos <p% + sin (p f) ± sin Ok z, 



(B.l) 
(B.2) 
(B.3) 



where the upper and lower signs in + cos Ok are for the inci- 
dent and reflected waves, respectively. The coordinates in which 
Eq. © is written are (x',y',z') (Fig. [TJ, defined according to re- 
lations y ' = B x k/\B x k\ and x' = y' x k. 

The electric field of the incoming ray with unit amplitude 
and polarization M' (M' =X or M' -O) can be written as 



(i) (0 (i) , (i) (i) 

'M' ~ c M'l e l T c M'2 e 2 ■ 



where c% . = e% 
reflected field is 



'XM' e X T 'OM ,t! Q 



„0) 



(B.4) 

According to Eqs. ( 1A.U and (IB.4t . the 



2 2 

V V 



m- 1 j= 1 



(B.5) 



The amplitudes r' XM , and r' QM , of the reflected-field components 
in the X- and O-modes, respectively, are given by the solution of 
the linear system 



(B.6) 



where c^. = e^ ■ ej 1 . Since the incident X- and O-modes are 
incoherent, the normal mode reflectivities in Eqs. d29l and ( l30l > 



L X1 L 01 


I r XM' 






( r & „® \ 

L X1 L 01 


M Jf) 

V C X2 c 02 > 


\ r OM' i 




rn j 


Ji) „® 
^ C X2 c 02 / 



are 

RMM' = \ r MM'\ 2 - 

According to Eq. ©, 



(i,r) _ 1A M X M- e j + y 1,1" e j 
C Mj 



(i.r) 



(B.7) 



(B.8) 



1 + l^ r) l 2 



where = Km(u\j). The explicit expressions for the scalar 
products in Eq. dB.8t are 

= sinOg sin^/sinai jr , (B.9) 
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y'lr ' e i ~ (cos 9b sin 9k ± sin 9b cos 9k cos <p)/ sin a^, (B.IO) 
ir' e 2 = (+cos9Bsin9k - sin 9 b cos 9k cos ip) /sin a UI , (B.ll) 
y'r' e 2 ~ + sin 9 b sin (p (cos 2 9k~ sin 2 6k)/ sin a UI . (B.12) 

Caution should be used when employing approximations for 
r m j in Eq. ( 1B.6I I if one of the normal modes is almost completely 
reflected, that is, Rmm> ~ 1- Such a situation occurs for the X- 
mode at E c i < E < Eq, if both k and B are close to normal 
(see Paper I). In this case the fit error may exceed (1 - Rxx) 
and result in Rxx > 1, which is unphysical. In particular, the 
fitting formulae presented below may occasionally give Rxx a 
few percent above 1 at very small 9b and 9k- In such instances 
one should truncate the mode-specific reflectivities, recovered 
from the fit, so as to fulfill the general condition < Rmm' < 1. 



B.2. Approximations for r mj 

For calculating Rmw according to Sect. IB. ll we use an analytic 
model of the complex reflectivity coefficients r m j, which agrees 
with the approximations derived in Sect. l2.3l and roughly repro- 



duces the computed dependences of r m j = \r m j\ exp(i(/> m j) on E 
for many characteristic geometry settings. For the squared mod- 
uli we use the following expressions: 



ln 2 |' 



\r 2 ir = 



in Region I, 

f E (l-Jm)(l-Jc) (B.13) 

+/c^l(1 - sin 6b sin 2 <p)/2 in Region II, 

in Region I, 

/ £ (1 + 7bi-27 b )(1-/c) (B.14) 

+JcRl sin 9k (1 - cos a)/2 in Region II, 



1-/1 



V22X 



1-/2 



l>-21 1 



(B.15) 



and/g = E/(E + Eq/2). The functions J\, J2, Jb, Jc, Jm, and/?£ 
are defined in Sect. 12.31 In the case of the free-ions model, small 
accidental discontinuities at the boundary of Region I are elimi- 
nated by truncating \r\ 1 1 2 and |r22| 2 from above by their values at 
E = E c i, Our approximations for the complex phases are 



n in Region I, 
nf L , if E > E c , 

E - E c - 
n + n — — otherwise, 



»22 = 



Eq - E c 

if E > E c , 
otherwise, 



nf L , 
<Pu 

-tt/2 in Region I, 

tt/2 - Infu if E > E c , 

E - (E ci + E c )/2 . 
n z otherwise, 



021 = 012 + 

where 



(B.16) 
(B.17) 

(B.18) 
(B.19) 



1 + exp 5 



E L -E 
El - Ec 
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